#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "clapack.h"
#include "richalgo.h"
#include "rhomx.h"
#include "spline.h"
#include "regularisation.h"
#include "materialfunction.h"
#include <plplot/plplot.h>
#include <time.h>



int main(int argc, char **argv) {
    int i;



    /**************************************/
    /* definition of timestep/size and gridsize/step */
    double dt         = 1;
    int    tmax       = 330;
    int    iter_t     = (int) tmax/dt;

    double dz         = 0.01;
    int    zmax       = 3;
    int iter_z        = (int) zmax/dz+1;

    /***************************************/
    /* definition of the material parameters*/
    double p0             = 325;                      /*capillary pressure startvalue*/
    double phi[3]         = {0.9 , 0.8889 , 0.95};		/*medium porosity*/

    double k0[3]          = {1500 , 90 , 12};        	/*permeability*/
    double material[3][3] = {{52.014 ,  0.715 ,   6.9},
                             {8.5 ,     0.96  ,   3.77},
                             {2.8 ,     1.674  ,  1.621}};
    double H0[iter_z];
    for(i = 0;i<iter_z;i++){
        if (i == 0){
            H0[i] = 0;
        }else{
        H0[i]        = -p0;
        }
    }
    double H01[iter_z];
    double H02[iter_z];
    double H004[iter_z];
    double H002[iter_z];
    memcpy(H01,H0,iter_z*sizeof(double));
    memcpy(H02,H0,iter_z*sizeof(double));
    memcpy(H002,H0,iter_z*sizeof(double));
    memcpy(H004,H0,iter_z*sizeof(double));

    /*********************************************/
    /*define rhomatrix*/
    double rhomatrix[iter_z][3];

    rhomxmat(rhomatrix, iter_z , zmax);
    /********************************************/



    /*************************************************/
    double rhomax01[iter_z][3];
    double rhomax02[iter_z][3];
    double rhomax002[iter_z][3];
    double rhomax004[iter_z][3];
    memcpy(rhomax01,rhomatrix,iter_z*sizeof(double)*3);
    memcpy(rhomax02,rhomatrix,iter_z*sizeof(double)*3);
    memcpy(rhomax002,rhomatrix,iter_z*sizeof(double)*3);
    memcpy(rhomax004,rhomatrix,iter_z*sizeof(double)*3);

    /***REGULARISATION******/
    double cornergap004 = 0.04;
    regular(rhomax004,cornergap004,dz,iter_z);
    double cornergap002 = 0.02;
    regular(rhomax002,cornergap002,dz,iter_z);
    double cornergap01 = 0.1;
    regular(rhomax01,cornergap01,dz,iter_z);
    double cornergap02 = 0.2;
    regular(rhomax02,cornergap02,dz,iter_z);

    /*************************************************/
    /*************************************************/
    /*******RICHARDALGORITHM**************************/
    clock_t prgstart,prgend;
    int hpicar;
    int hpicar01;
    int hpicar02;
    int hpicar002;
    int hpicar004;


    prgstart = clock();
    hpicar = richardsalgo(H0,1,iter_z,dz,dt,material,k0,phi,rhomatrix,iter_t);
    prgend = clock();
    for(i = 0;i<iter_z;i++){
        printf("H %f \n",H0[i]);
    }
//    double time0 = (double)(prgend-prgstart) / CLOCKS_PER_SEC;
//    printf("time %f \n",time0);
//    prgstart = clock();
//    hpicar01 = richardsalgo(H01,1,iter_z,dz,dt,material,k0,phi,rhomax01,iter_t);
//    prgend = clock();
//    double time01 = (double)(prgend-prgstart) / CLOCKS_PER_SEC;

//    prgstart = clock();
//    hpicar02 = richardsalgo(H02,1,iter_z,dz,dt,material,k0,phi,rhomax02,iter_t);
//    prgend = clock();
//    double time02 = (double)(prgend-prgstart) / CLOCKS_PER_SEC;

//    prgstart = clock();
//    hpicar002 = richardsalgo(H002,1,iter_z,dz,dt,material,k0,phi,rhomax002,iter_t);
//    prgend = clock();
//    double time002 = (double)(prgend-prgstart) / CLOCKS_PER_SEC;

//    prgstart = clock();
//    hpicar004 = richardsalgo(H004,1,iter_z,dz,dt,material,k0,phi,rhomax004,iter_t);
//    prgend = clock();
//    double time004 = (double)(prgend-prgstart) / CLOCKS_PER_SEC;

//    double cputime[5] = {time0,time002,time004,time01,time02};
//    double cornergaparray[5] = {0.0 , cornergap002 ,cornergap004,cornergap01,cornergap02};
//    double picariter[5] = {hpicar, hpicar002, hpicar004, hpicar01,hpicar02};


//    double izvec[iter_z];
//    for(i = 0;i<iter_z;i++){
//        izvec[i] = i*dz;
//    }

//    PLFLT xmin,xmax,ymin,ymax;//x[5] = {0,1,2,3,4},y[5]={0,1,4,9,16},z[5] = {1,2,1,2,1};
//    PLINT just;
//    PLINT axis;
//    PLINT n=iter_z;

//    PLINT color2 = 2;
//    PLINT color3 = 2;
//    PLINT color4 = 4;
//    PLINT color5 = 5;
//    PLINT color9 = 9;

//    //PLINT style = 9;




//    just = 0;
//    axis = 0;
//    xmin = 0.0;
//    xmax = 4;
//    ymin = -p0;
//    ymax = 0.0 ;



//    c_plsdev("psc"); // colored postscript
//    c_plsfnam("pressurehead.ps");
//    c_plinit();  // start plplot object
//    c_plenv(xmin, xmax, ymin, ymax, just, axis );

//    c_pllab( "(depth)", "(pressure head)", "Pressurehead dz = 0.001 dt = 0.4");

//    c_plcol0(color2);
//    c_plline( n, izvec, H0);
//    c_plcol0(color3);
//    c_plline( n, izvec, H02);
//    c_plcol0(color4);
//    c_plline( n, izvec, H01);
//    c_plcol0(color5);
//    c_plline( n, izvec, H002);
//    c_plcol0(color9);
//    c_plline( n, izvec, H004);

//    c_plend();
//    just = 0;
//    axis = 0;
//    xmin = 0.0;
//    xmax = 0.5;
//    ymin = 0.0;
//    ymax = cputime[0]+2 ;
//    PLINT m = 5;
//    PLINT style = 9;

//    c_plsdev("psc"); // colored postscript
//    c_plsfnam("cputime.ps");
//    c_plinit();  // start plplot object
//    c_plenv(xmin, xmax, ymin, ymax, just, axis );

//    c_pllab( "(cornergap)", "(time)", "Comparison Time");
//    c_plcol0(color9);
//    c_plpoin(m,cornergaparray,cputime,style);
//    c_plend();

//    just = 0;
//    axis = 0;
//    xmin = 0.0;
//    xmax = 0.5;
//    ymin = 0.0;
//    ymax = picariter[0] ;
//    PLINT k = 5;


//    c_plsdev("psc"); // colored postscript
//    c_plsfnam("picariter.ps");
//    c_plinit();  // start plplot object
//    c_plenv(xmin, xmax, ymin, ymax, just, axis );

//    c_pllab( "(cornergap)", "(Iterations Picar)", "Comparison Iterations");
//    c_plcol0(color9);
//    c_plpoin(k,cornergaparray,picariter,style);
//    c_plend();



    return 0;
}//end main


